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We present molecular dynamics simulations of aggregation kinetics in a colloidal suspension mod- 
eled as a highly asymmetric binary mixture. Starting from a configuration with largely uncorrelated 
colloidal particles the system relaxes by coagulation-fragmentation dynamics to a structured state 
of low-dimensionality clusters with an exponential size distribution. The results show that short- 
range repulsive interactions alone can give rise to so-called cluster phases. For the present model 
and probably other, more common colloids, the observed clusters appear to be equilibrium phase 
fluctuations induced by the entropic inter-colloidal attractions. 
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Colloidal suspensions are common in biological and geological systems and their unusual properties are intensely 
exploited in many industrial applications, from pharmaceutics to food processing. Although some, e.g. noble metal 
colloids, have been produced and used for centuries, the systematic study of these systems is a more recent endeavor 
. One of the remarkable features of colloids is the possibility of tunin g th e colloidal interactions . This multiplies 
the already rich physical behavior exhibited by colloids |^lj>IMEHBlEllflllllll3>ll3] an d significantly increases their 
technological potential, but it also underscores the need for understanding the general principles and mechanisms that 
govern their equilibrium and kinetic behavior. The low volume fractions kinetic arrest leading to gelation in weakly 
attractive hard-sphere colloids has been found for example to share many of the features of the high volume fractions 
colloidal glass transition Q . The aggregation of colloids and the ensuing highly structured states such as gels are 
a direct consequence of the short-ranged attractions typically encountered in these systems, known as depletion or 
entropic interactions The role of repulsions however is also important, and recently much attention has been 

devoted to elucidating the effect of long-range repulsions in the formation of cluster phases and the gelation transition 
0, flfi. UtI ] . Understanding clustering effects and the mechanisms that control them is important conceptually as well 
as practically in systems ranging from biological to engineered ones, e.g. nanofluids. In the present paper we explore 
using molecular dynamics (MD) simulations the characteristics of aggregation kinetics in a model colloidal suspension 
with short-range repulsions that explicitly includes both colloidal and solvent particles. We find that the system 
evolution is well described as coagulation-fragmentation dynamics surprisingly similar to a simple polymerization- 
depolymerization process. The final state is a structured phase of low-dimensionality clusters with an exponential 
size distribution. These properties appear to be all due to the bond-like character of the induced short-ranged 
inter-colloidal attractions. 

The simplest instance of a colloidal suspension is perhaps a binary hard-sphere mixture with significant size asym- 
metry. Dijkstra et al. 18] have shown that the phase behavior of such a system can be accurately predicted by 
integrating out the effect of the small particles and introducing an effective depletion potential between the large 
(colloidal) particles. This approach provides significant insight into the structure of the phase diagram and simplifies 
considerably the computational study of the mixture. However, it may not be necessarily appropriate for the study 
of non-equilibrium processes in colloids. The motion of colloidal particles suspended in viscous liquids is generally 
coupled, even at large separations, by the induced suspending-liquid flows. These so-called hydrodynamic interactions 
are believed for example to play an important role in the aggregation kinetics . In the following we consider and 
fully model using MD simulations a system related to the binary hard-sphere mixture, consisting of two types of 
particles, 1 - colloid and 2 - solvent, with equal masses to. The interactions between the particles are based on the 
inverse-12, 'soft-sphere' potential, whose properties have been well studied [2(| 

u(r)=e^) 12 (1) 

, and which we truncate and shift at r/a = 2. (We also define u(r) — oo for r < 0.) The interactions are: 

wu(r) =u(r-2R c ) (2a) 
uia(r) = u(r - (2b) 
u 2 2{r) = u(r) (2c) 

Similar potentials, that take into account the 'size' of the colloidal particles by introducing a hard core radius R c , have 
been employed before to model suspensions |2lj . For temperatures fc^T ~ e the effective diameters corresponding 
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to the above interactions should be well approximated by er 2 = a, a 12 — Rc + c, and o\ = a c = 2R C + a, and 
satisfy additivity, a\2 — {cri + 02)/2. The particular system that we focus on has a colloid-solvent 'diameter' ratio 
a c /a = 5, colloid 'volume fraction' <p c — irn c a c 3 /6 — 0.1 and solvent 'volume fraction' <f> = nna 3 /6 — 0.37; n c and 
n are the corresponding number densities, n c = N c /V, n — N/V. We perform MD simulations of this system in 
the microcanonical (NVE) ensemble with iV c = 200, N ~ 10 5 and an av erag e temperature set to fcgT = e. Such 
simulations of binary mixtures with large size ratios have been done before |22(, but not, to our knowledge, for colloid 
volume fractions as small as the present one. 

The simulations are initiated from a configuration of non-overlapping (as defined by the effective diameters), but 
otherwise randomly distributed colloid and solvent particles, and the system is equilibrated at fc^T = e using a 
combination of velocity rescaling and Anderson thermostatting. Thermal equilibrium is achieved over time scales 
much shorter than the scales over which the colloidal particles move over any significant distances, about 5000 MD 
steps, leading to a state corresponding to a fully equilibrated solvent and essentially uncorrelated colloidal particles. 
A largely similar configuration would be expected for example in the case of a hard-sphere colloid shortly after a rapid 
quench from a high temperature state. The system evolution is subsequently followed using conventional (constant 
energy) MD for over 5 x 10 5 time steps. 

The evolving configuration of the large (colloidal) particles during the simulations is best characterized by using a 
standard aggregation criterion, that places two particles in the same cluster when their relative distance is smaller 
than a certain value r^. For the present work we set r& = 1.05ct c , but our conclusions are independent of this value 
in a fairly wide range. The application of this criterion yields a time dependent cluster size distribution Nk(t), or in 
scaled form nk(t) = Nk(t)/N c , where N c — J^k kNk{t) is the (constant) total number of particles. Wc show in Fig. 1 
the evolution of the (normalized) number of monomers ni(t), along with the inverse of the (normalized) total number 
of clusters, si(t) = '}2 lk kNk{t)/'^2 lk Nk(t)- 1 s\ can also be interpreted as a measure of the average cluster size. The 
number of monomers in the system steadily decreases, while the cluster size increases, with both quantities apparently 
reaching at late times well defined plateaus. Visual inspection of the clustering dynamics reveals tenuous structures 
reminiscent of branched polymers - see Fig. 2 - whose evolution is driven by aggregation and breakup processes. 
These features suggest that general concepts developed to describe coagulation-fragmentation kinetics [3 E] mav 
also be relevant to the present system. The starting point for this analysis is the generalized Smoluchowski equation 
for the cluster size distribution: 

drik 1 00 

— = - [KijUiUj - FijUi+j] - ^yKkjUurij - Fkjiik+j] (3) 

i+j=k j = l 

with coagulation and breakup rate coefficients {-fQj} and {-Fy}, i,j — 1, 00. The cluster population rik(t) whose 
evolution is described by the above equations is generally assumed to satisfy a dynamical scaling relation |2.'il Eij l 

n k {t) = s{t)- 2 cp(x) (4) 

where x = k/s(i), s(t) is a characteristic cluster size and (j)(x) a scaling function dependent on the details of process, 
i.e. the coefficients {Kij} and {F^}. The applicability of Eq. 0]to systems that reach equilibrium, i.e. where s(t) goes 
to a constant value when t —> 00, appears to be somewhat more limited |25j than initially assumed. Nevertheless, 
under certain reasonable assumptions for the coagulation and breakup rates and in particular for a classical model of 
polymerization-depolymerization with constant coefficients |26| dynamical scaling holds exactly. For small deviations 
from equilibrium it can also be shown [241 that dynamical scaling implies that s(t), and therefore si(t), relaxes 
exponentially to its final, equilibrium value |24j. Furthermore, for the case of Ref. |26( the product n\{t)s\(t) is time 
independent, ni(t)s\(t) = 1, while in general should increase only slowly in a narrow range. 

For the asymmetric binary mixture considered here as an archetype of a colloidal suspension we find that s±(t) 
relaxes to a very good approximation in an exponential fashion, and moreover that the relation ni(t)s\(t) = 1 is 
very close to being satisfied - see Fig. 1. This makes the behavior of this simple model, containing only short-range 
repulsive potentials, surprisingly similar to that of common associating systems. As a test of dynamical scaling we 
consider the usual measure of the mean cluster size H3,|3; s (t) = Ylk k 2 nk(t), which depends linearly on si(i) for 
the model of Blatz and Tobolsky j2|| , and in general should be proportional with it when s is large enough. This is 
also rather well satisfied, see Fig. 1 (inset). 

The morphology of the typical clusters - see Fig. 2 - closely resembles the open structures usually encountered in 
colloidal aggregation 1 13. [271 ■ Although the observed assemblies are not nearly as large as those obtained upon 
irreversible flocculation |27|. their dimensionality can be well quantified by considering their radius of gyration 
kRg(k) = (J2i=ii( x i ~ x cm) 2 + (Vi ~ Ucm) 2 + ( z i — z cm) 2 }), where i denotes the particles in a cluster of size k, 
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CM is the cluster center of mass, and the average is taken over the final population of clusters of size k. The resulting 
radii - see Fig. 3 - obey R g (k) oc k 1 ' D f, with Df ~ 2. The fractal dimension Df is slightly higher than the one 
observed for typical colloidal aggregates [27| , but in agreement with that of randomly branched polymers [28| . 

The cluster size distribution in the final state - n^Jk), see Fig. 4 - compares well with the distributions obtained 
in model coagulation-fragmentation dynamics [23I |26| and may be reasonably described at large cluster sizes as a 
decaying exponential. From the point of view of the kinetics this appears to be a consequence of an approximately 
cluster-independent break-up rate |23l |26| , which is consistent with the low cluster dimensionality and the equivalence 
for this system of the 'bonds' that hold together structures of different sizes. In fact, thermodynamic arguments 
applicable to equilibrium polymers systems 1291 yield a similar prediction, while general concepts pertaining to phase 
fluctuations (and nucleation) in solutions [3J|, would also predict an exponentially decreasing abundance of large 
clusters for single-phase equilibrium close to a two-phase coexistence region. Nevertheless, since the size of our system 
in terms of the number of large particles modeled is fairly small, we further address this issue by considering an 
equivalent one-component model with effective interactions. 

To this end we consider the structure of the mixture, in particular the colloid-colloid pair correlation function, 
9n( r ) - Fig. 5. This function exhibits the sharp nearest-neighbor peak also encountered in the effective one-component 
modeling of asymmetric binary hard-sphere mixtures |lSj. which is due to the strong, entropically induced attractions. 
However, gu(r) is also considerably more structured than such one-component correlation functions due to the explicit 
presence of the solvent. For our own single component modeling we consider particles interacting through a potential 
V e ff(r) defined by (3V e ff(r) = — ln[gn(r)], i.e. the potential of mean force, but which is truncated after its first 
maximum at a distance r rn satisfying V e ff(r m ) = 0. We expect the resulting interaction to be an appropriate effective 
pair potential for the rather low colloid volume fraction system studied here. V e /f has a fairly deep well followed 
by a barrier - see Fig. 5 (inset) - and is analogous with the effective potentials used to study the phase diagram of 
asymmetric binary hard-sphere mixtures jig ■ MD simulations with 864 particles interacting through this potential at 
a temperature fcsT = e were carried out. The system requires very long equilibration times, of the order of 10 6 steps, 
which we then follow by standard (NVE) MD to determine the fluid structure. The main features of the resulting 
pair correlation function - Fig. 5 - are not surprising. The sharp nearest neighbor peak of the binary mixture along 
with the next (colloid-separated) nearest neighbor maximum are very well reproduced, but the intermediate structure 
due to the small particles is missing. To further elucidate the one-component system structure and compare it with 
that of the binary mixture we perform the same clustering analysis as before, using the same distance r&. This yields 
the cluster population shown in Fig. 4, which agrees well with the fully modeled mixture result and exhibits a rather 
clear exponential decay at large cluster sizes. Somewhat surprisingly, the dimensionality of the clusters so obtained 
is also in very good agreement with that of the clusters encountered in the two-component simulation of the mixture. 
This appears to cast some doubt on the role played by hydrodynamic interactions in the formation of open structures 

The present MD results and analysis suggest that in strongly asymmetric binary mixtures exclusion-like interactions 
alone can give rise to highly structured, so-called cluster phases. For the binary mixture studied here as a paradigm 
of a colloidal suspension the clustering kinetics is dominated by two-body coagulation and fragmentation processes. 
This may be expected to hold with increasing accuracy as the colloidal volume fraction decreases. The resulting 
colloidal clusters are tenuous, open structures that can be classified as randomly branched polymers, which is largely 
a consequence of the very short range of the induced inter-colloidal attractions and not hydrodynamic interactions. In 
the final state the clusters are exponentially distributed in size for both the fully modeled mixture and an equivalent 
effective one-component fluid. This is consistent with single-phase equilibrium of the mixture close to a two-phase 
region of its phase diagram. The observed cluster 'phase' is then simply a signature of equilibrium phase fluctuations, 
controllable for example by temperature and pressure changes. This may also be the case with other, more standard 
examples of colloids where clusters form through reversible aggregation and are quite small in size [j] ■ When considered 
together with the observed properties of the aggregation kinetics these features lend support to the idea [24| that the 
relaxation dynamics in such systems is well described as coagulation- fragmentation kinetics of phase fluctuations. 

This work was performed under the auspices of the U. S. Department of Energy by University of California Lawrence 
Livermore National Laboratory under Contract No. W-7405-Eng-48. 
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FIG. 1: Time evolution of the (normalized) number of monomers - ni (bottom dotted line), inverse of the (normalized) total 
number of clusters - si (top dotted line), and nis? (middle dotted line); solid line is an exponential relaxation fit; to = o(m/ e) 1 / 2 . 
Inset: characteristic cluster size s (see text) versus si. 




FIG. 2: A cluster of 33 colloidal particles. 
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FIG. 3: Radius of gyration R g (see text) as a function of cluster size k: from simulations (dashed line) and power law 
(R g (x k 1/D i) fit (solid line); D f = 2. 




FIG. 4: Cluster size distribution from MD simulations of the binary mixture (diamonds), simulations of the one-component 
fluid (dotted line) and exponential decay fit (solid line). 




FIG. 5: Colloid-colloid pair correlation function of the binary mixture (solid line) and effective one-component system (dotted 
line). Inset: interaction potential V e ff (see text) of the one-component system; /3 = 1/fcsT. 



